import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
import healpy as hp
import pixell as px
from pixell import reproject
from classy import Class
pi = np.pi
plt.style.use(['mysize_meeting'])
fig = plt.figure()
fw = fig.get_figwidth()
fh = fig.get_figheight()
<Figure size 720x540 with 0 Axes>
cmap = plt.cm.RdBu
c__256 = 'C0'
c__512 = 'C1'
c_1024 = 'C2'
c_proj = 'C2'
c_fgen = 'C3'
The resolution of a full-sky map is governed by the parameter nside.
This is also where I specify the image size used throughout this notebook.
res2048_arcmin = hp.nside2resol(2048, arcmin=True)
res1024_arcmin = hp.nside2resol(1024, arcmin=True)
res512_arcmin = hp.nside2resol(512, arcmin=True)
res256_arcmin = hp.nside2resol(256, arcmin=True)
print(res2048_arcmin)
print(res1024_arcmin)
print(res512_arcmin)
print(res256_arcmin)
1.717743205908703 3.435486411817406 6.870972823634812 13.741945647269624
imgsize_deg = 100 # in degrees
imgsize_arcmin = imgsize_deg * 60 # in arcmin
imgsize_rad = imgsize_deg / 180 * pi # in radians
pixsize_256 = int(np.ceil(imgsize_arcmin / res256_arcmin))
pixsize_512 = int(np.ceil(imgsize_arcmin / res512_arcmin))
pixsize_1024 = int(np.ceil(imgsize_arcmin / res1024_arcmin))
imgsize_arcmin = pixsize_1024 * res1024_arcmin
print("img: %.2f deg = %.f arcmin" % (imgsize_arcmin/60, imgsize_arcmin))
print("pix: %d pixels" % pixsize_1024)
print("res: %.3f deg = %.2f arcmin" % (res1024_arcmin/60, res1024_arcmin))
img: 100.03 deg = 6002 arcmin pix: 1747 pixels res: 0.057 deg = 3.44 arcmin
try:
cosmo.struct_cleanup()
except:
pass
cosmo = Class()
cosmo.set({
'output': 'tCl lCl',
'modes': 's',
'omega_b': 0.022383,
'omega_cdm': 0.12011,
'100*theta_s': 1.041085,
'tau_reio': 0.0543,
'A_s': np.exp(3.0448)*1e-10,
'n_s': 0.96605,
'm_ncdm': 0.06,
'N_ncdm': 1,
'N_ur': 2.0328,
'lensing': 'yes',
'l_max_scalars': 20000
})
True
%%time
cosmo.compute()
CPU times: user 1min 30s, sys: 3.28 s, total: 1min 33s Wall time: 7.33 s
ell = cosmo.lensed_cl()['ell']
lmax = ell[-1]
norm = ell * (ell+1) / (2*np.pi)
cl_cmb = cosmo.lensed_cl()['tt'] * cosmo.T_cmb()**2 / 1e-6**2
cl_pwl = ell**2 / (1e9 + ell**4) * 2e8
cl_pwl[1:] /= norm[1:]
def power_axes(ax1, ax2):
ax1.set_xlabel("$\ell$")
ax2.set_xlabel("$\ell$")
ax1.set_ylabel("$C_\ell\ \ell (\ell+1) / 2\pi\ /\ \mu\mathrm{K^2}$")
ax2.set_ylabel("$C_\ell\ \ell (\ell+1) / 2\pi\ /\ \mu\mathrm{K^2}$")
ax1.set_xlim(-100, 3000)
ax2.set_xlim(1, 10000)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_cmb*norm, c='k', label="CMB power spectrum (from CLASS with P18 best-fit parameters)")
ax1.plot(ell, cl_pwl*norm, c='0.5', label="custom power-laws power spectrum")
ax1.legend(handlelength=0, labelcolor='linecolor')
ax2 = plt.subplot(212)
ax2.loglog(ell, cl_cmb*norm, c='k' )
ax2.loglog(ell, cl_pwl*norm, c='0.5')
ax2.set_ylim((cl_pwl*norm)[2]/10, (cl_pwl*norm).max()*4)
power_axes(ax1, ax2)
%%time
alm_cmb = hp.synalm(cl_cmb, lmax=lmax)
alm_pwl = hp.synalm(cl_pwl, lmax=lmax)
CPU times: user 30.5 s, sys: 1.8 s, total: 32.3 s Wall time: 33.8 s
%%time
map_cmb_256 = hp.alm2map(alm_cmb, nside=256, lmax=lmax)
map_cmb_512 = hp.alm2map(alm_cmb, nside=512, lmax=lmax)
map_cmb_1024 = hp.alm2map(alm_cmb, nside=1024, lmax=lmax)
map_pwl_256 = hp.alm2map(alm_pwl, nside=256, lmax=lmax)
map_pwl_512 = hp.alm2map(alm_pwl, nside=512, lmax=lmax)
map_pwl_1024 = hp.alm2map(alm_pwl, nside=1024, lmax=lmax)
CPU times: user 15min 14s, sys: 705 ms, total: 15min 14s Wall time: 1min 1s
vmax = max(abs(map_cmb_1024.min()), abs(map_cmb_1024.max()))
hp.mollview(map_cmb_1024, min=-vmax, max=vmax, cmap=cmap, title="CMB map")
vmax = max(abs(map_pwl_1024.min()), abs(map_pwl_1024.max()))
hp.mollview(map_pwl_1024, min=-vmax, max=vmax, cmap=cmap, title="custom power-law map")
%%time
cl_cmb_out_256 = hp.anafast(map_cmb_256, lmax=lmax, use_pixel_weights=True)
cl_cmb_out_512 = hp.anafast(map_cmb_512, lmax=lmax, use_pixel_weights=True)
cl_cmb_out_1024 = hp.anafast(map_cmb_1024, lmax=lmax, use_pixel_weights=True)
cl_pwl_out_256 = hp.anafast(map_pwl_256, lmax=lmax, use_pixel_weights=True)
cl_pwl_out_512 = hp.anafast(map_pwl_512, lmax=lmax, use_pixel_weights=True)
cl_pwl_out_1024 = hp.anafast(map_pwl_1024, lmax=lmax, use_pixel_weights=True)
CPU times: user 15min 7s, sys: 3.72 s, total: 15min 11s Wall time: 1min 13s
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_cmb_out_256 * norm, c=c__256, label="from map with nside=256")
ax1.plot(ell, cl_cmb_out_512 * norm, c=c__512, label="from map with nside=512")
ax1.plot(ell, cl_cmb_out_1024 * norm, c=c_1024, label="from map with nside=1024")
ax1.plot(ell, cl_cmb * norm, c='k', label="CMB input Cl")
ax1.set_ylim(0, max(cl_cmb*norm)*1.2)
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor')
ax2 = plt.subplot(212)
ax2.loglog(ell, cl_cmb_out_256 * norm, c=c__256)
ax2.loglog(ell, cl_cmb_out_512 * norm, c=c__512)
ax2.loglog(ell, cl_cmb_out_1024 * norm, c=c_1024)
ax2.loglog(ell, cl_cmb * norm, c='k')
ax2.set_ylim((cl_cmb*norm)[2]*1e-2, max(cl_cmb*norm)*4)
power_axes(ax1, ax2)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(ell, cl_pwl_out_256 * norm, c=c__256, label="from map with nside=256")
ax1.plot(ell, cl_pwl_out_512 * norm, c=c__512, label="from map with nside=512")
ax1.plot(ell, cl_pwl_out_1024 * norm, c=c_1024, label="from map with nside=1024")
ax1.plot(ell, cl_pwl * norm, c='k', label="custom power-law input Cl")
ax1.set_ylim(0, max(cl_pwl*norm)*1.2)
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor')
ax2 = plt.subplot(212)
ax2.loglog(ell, cl_pwl_out_256 * norm, c=c__256)
ax2.loglog(ell, cl_pwl_out_512 * norm, c=c__512)
ax2.loglog(ell, cl_pwl_out_1024 * norm, c=c_1024)
ax2.loglog(ell, cl_pwl * norm, c='k')
ax2.set_ylim((cl_pwl*norm)[2]*1e-2, max(cl_pwl*norm)*4)
power_axes(ax1, ax2)
def fourier_power_spectrum(img_r, img_len, bin_k):
"""
Compute the Fourier power spectrum given a 2-dimensional real-space map.
Parameters
----------
img_r: arraylike
2-dimensional real-space map
img_len: float
image size
(units define output units of wavenumber and power spectrum)
bin_k: int
number of bins in Fourier space
Returns
-------
k: np.ndarray
wavenumber in units of inverse [img_len]
pk: np.ndarray
power spectrum in units of [img_len]**2
"""
assert np.ndim(img_r) == 2, 'img_r is not 2D'
img_pix = np.size(img_r, axis=0) # pixel number per axis
img_dim = np.ndim(img_r) # image dimension
# This first 'paragraph' is to create masks of indices corresponding to
# one Fourier bin each.
freq = np.fft.fftfreq(n=img_pix, d=img_len/img_pix) * 2*np.pi
rfreq = np.fft.rfftfreq(n=img_pix, d=img_len/img_pix) * 2*np.pi
kx, ky = np.meshgrid(freq, rfreq, indexing='ij')
k_abs = np.sqrt(kx**2. + ky**2.)
# The following complicated line is actually only creating a 1D array
# spanning k-space logarithmically from minimum k_abs to maximum. To
# start slightly below the minimum and finish slightly above the maximum
# I use ceil and floor. To ceil and floor not to the next integer but to
# the next 15th digit, I multiply by 1e15 before flooring and divide by
# the same amount afterwards. Since the ceiled/floored value is actually
# the exponent used for the logspace, going to the next integer would be
# too much, which is why I do this ceil/floor to digit procedure.
k_log = np.logspace(np.floor(np.log10(np.min(k_abs[1:]))*1.e15)/1.e15, np.ceil(np.log10(np.max(k_abs[1:]))*1.e15)/1.e15, bin_k)
img_k = np.fft.rfftn(np.fft.fftshift(img_r)) * (img_len/img_pix)**img_dim
pk = np.empty(np.size(k_log)-1)
for i in range(np.size(k_log)-1):
mask = (k_abs >= k_log[i]) & (k_abs < k_log[i+1])
pk[i] = np.mean(np.abs(img_k[mask])**2.) / img_len**img_dim
k = (k_log[1:] + k_log[:-1]) / 2
return k, pk
Using the pixell package as mentioned in the ACT example notebooks to project the full-sky map to a flat cut-out.
from pixell import enmap, utils
box = np.array([[-imgsize_arcmin//2, imgsize_arcmin//2], [imgsize_arcmin//2, -imgsize_arcmin//2]]) * utils.arcmin
shape256, wcs256 = enmap.geometry(pos=box, res=res256_arcmin*utils.arcmin, proj='car')
shape512, wcs512 = enmap.geometry(pos=box, res=res512_arcmin*utils.arcmin, proj='car')
shape1024, wcs1024 = enmap.geometry(pos=box, res=res1024_arcmin*utils.arcmin, proj='car')
imap256 = enmap.zeros((3,) + shape256, wcs=wcs256)
imap512 = enmap.zeros((3,) + shape512, wcs=wcs512)
imap1024 = enmap.zeros((3,) + shape1024, wcs=wcs1024)
imap1024.shape
(3, 1747, 1747)
imap1024.wcs
car:{cdelt:[-0.05726,0.05726],crval:[-0.008333,0],crpix:[874.1,874]}
%%time
fmap_cmb_256 = reproject.enmap_from_healpix(map_cmb_256, shape=(pixsize_256, pixsize_256), wcs=imap256.wcs)
fmap_cmb_512 = reproject.enmap_from_healpix(map_cmb_512, shape=(pixsize_512, pixsize_512), wcs=imap512.wcs)
fmap_cmb_1024 = reproject.enmap_from_healpix(map_cmb_1024, shape=(pixsize_1024, pixsize_1024), wcs=imap1024.wcs)
fmap_pwl_256 = reproject.enmap_from_healpix(map_pwl_256, shape=(pixsize_256, pixsize_256), wcs=imap256.wcs)
fmap_pwl_512 = reproject.enmap_from_healpix(map_pwl_512, shape=(pixsize_512, pixsize_512), wcs=imap512.wcs)
fmap_pwl_1024 = reproject.enmap_from_healpix(map_pwl_1024, shape=(pixsize_1024, pixsize_1024), wcs=imap1024.wcs)
Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting Preparing SHT T -> alm float64 complex128 Computing pixel positions Computing rotated positions Projecting CPU times: user 57.5 s, sys: 980 ms, total: 58.5 s Wall time: 13.7 s
vmax = max(abs(fmap_cmb_1024[0].min()), abs(fmap_cmb_1024[0].max()))
plt.title("flat projection from CMB map\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_cmb_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7fa20efbaa90>
vmax = max(abs(fmap_pwl_1024[0].min()), abs(fmap_pwl_1024[0].max()))
plt.title("flat projection from custom power-law map\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_pwl_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7fa20f0c1040>
fmap_cmb_256_k, fmap_cmb_256_pk = fourier_power_spectrum(fmap_cmb_256[0], img_len=imgsize_arcmin, bin_k=256)
fmap_cmb_512_k, fmap_cmb_512_pk = fourier_power_spectrum(fmap_cmb_512[0], img_len=imgsize_arcmin, bin_k=256)
fmap_cmb_1024_k, fmap_cmb_1024_pk = fourier_power_spectrum(fmap_cmb_1024[0], img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_256_k, fmap_pwl_256_pk = fourier_power_spectrum(fmap_pwl_256[0], img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_512_k, fmap_pwl_512_pk = fourier_power_spectrum(fmap_pwl_512[0], img_len=imgsize_arcmin, bin_k=256)
fmap_pwl_1024_k, fmap_pwl_1024_pk = fourier_power_spectrum(fmap_pwl_1024[0], img_len=imgsize_arcmin, bin_k=256)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount) /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount)
def power_axes_k(ax1, ax2):
ax1.set_xlabel("$k$")
ax2.set_xlabel("$k$")
ax1.set_ylabel("$P(k)\ k^2 / (2\pi)^2$")
ax2.set_ylabel("$P(k)\ k^2 / (2\pi)^2$")
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(fmap_cmb_512_k, fmap_cmb_512_k**2/(2*pi)**2 * fmap_cmb_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(fmap_cmb_256_k, fmap_cmb_256_k**2/(2*pi)**2 * fmap_cmb_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the $\\bf{projection}$ to a flat map \nof the CMB power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(fmap_cmb_512_k, fmap_cmb_512_k**2/(2*pi)**2 * fmap_cmb_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(fmap_cmb_256_k, fmap_cmb_256_k**2/(2*pi)**2 * fmap_cmb_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(fmap_pwl_512_k, fmap_pwl_512_k**2/(2*pi)**2 * fmap_pwl_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(fmap_pwl_256_k, fmap_pwl_256_k**2/(2*pi)**2 * fmap_pwl_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the $\\bf{projection}$ to a flat map \nof the custom power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(fmap_pwl_512_k, fmap_pwl_512_k**2/(2*pi)**2 * fmap_pwl_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(fmap_pwl_256_k, fmap_pwl_256_k**2/(2*pi)**2 * fmap_pwl_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
class GaussianFluctuationMap:
"""
Generator class for anisotropy maps from of Gaussian fluctuations for given angular power spectra.
"""
def __init__(self, ell, cl, len_rad, res_rad):
"""
Parameters
----------
ell: list or np.ndarray
multipole moment (needs to have the same shape as cl)
cl: list or np.ndarray
angular power spectrum (needs to have the same shape as ell)
len_rad: float
desired image size in radians
res_rad: float
desired image resolution in radians
"""
self.len_rad = len_rad # size of map in radians on the sky
self.res_rad = res_rad # resolution of the map in radians
self.len_pix = self.rad2pix(len_rad) # size of map in pixels
self.ell = ell[2:]
self.cl = cl[2:]
self.l2cl = interp1d(ell, cl)
def rad2pix(self, rad):
"""
Convert radians in sky to pixels in an image at a given resolution.
Parameters
----------
rad: float or np.ndarray
Returns
-------
pix: float or np.ndarray
number of pixels corresponding to the given angel `rad` at the given resolution
"""
pix = rad / self.res_rad
if type(pix) == float:
return int(pix)
else:
return pix.astype(int)
@staticmethod
def k2l_rad(k_irad):
"""
Convert wavenumber k to multipolemoment l.
Parameters
----------
k_irad: float or np.ndarray
wavenumber k in inverse radians
Returns
-------
multipolemoment `ell` (dimensionless)
"""
return 2 * pi * k_irad
def cl2map(self):
"""
Generate a random realization from the given specifications.
"""
k_rfreq = np.fft.rfftfreq(self.len_pix, d=self.res_rad)
k_freq = np.fft.fftfreq(self.len_pix, d=self.res_rad)
kx, ky = np.meshgrid(k_freq, k_rfreq, indexing='ij')
k_abs = np.sqrt(kx**2 + ky**2)
cl = self.l2cl(self.k2l_rad(k_abs))
# create the random map in k-space:
k_img = np.sqrt(cl * self.len_rad**2/2) * (np.random.randn(*cl.shape) +
1j * np.random.randn(*cl.shape))
# convert the map to real space (using an inverse Fourier transform):
r_img = np.fft.fftshift(np.fft.irfft2(k_img)) / self.res_rad**2
# devision by res to account for dx dy in FT, **2 for 2D
return r_img
gm_cmb_256 = GaussianFluctuationMap(ell=ell,
cl=cl_cmb,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res256_arcmin/60/180*pi)
gm_cmb_512 = GaussianFluctuationMap(ell=ell,
cl=cl_cmb,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res512_arcmin/60/180*pi)
gm_cmb_1024 = GaussianFluctuationMap(ell=ell,
cl=cl_cmb,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res1024_arcmin/60/180*pi)
gm_pwl_256 = GaussianFluctuationMap(ell=ell,
cl=cl_pwl,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res256_arcmin/60/180*pi)
gm_pwl_512 = GaussianFluctuationMap(ell=ell,
cl=cl_pwl,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res512_arcmin/60/180*pi)
gm_pwl_1024 = GaussianFluctuationMap(ell=ell,
cl=cl_pwl,
len_rad=imgsize_arcmin/60/180*pi,
res_rad=res1024_arcmin/60/180*pi)
gmap_cmb_256 = gm_cmb_256.cl2map()
gmap_cmb_512 = gm_cmb_512.cl2map()
gmap_cmb_1024 = gm_cmb_1024.cl2map()
gmap_pwl_256 = gm_pwl_256.cl2map()
gmap_pwl_512 = gm_pwl_512.cl2map()
gmap_pwl_1024 = gm_pwl_1024.cl2map()
gmap_cmb_256_k, gmap_cmb_256_pk = fourier_power_spectrum(gmap_cmb_256, img_len=imgsize_arcmin, bin_k=512)
gmap_cmb_512_k, gmap_cmb_512_pk = fourier_power_spectrum(gmap_cmb_512, img_len=imgsize_arcmin, bin_k=512)
gmap_cmb_1024_k, gmap_cmb_1024_pk = fourier_power_spectrum(gmap_cmb_1024, img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_256_k, gmap_pwl_256_pk = fourier_power_spectrum(gmap_pwl_256, img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_512_k, gmap_pwl_512_pk = fourier_power_spectrum(gmap_pwl_512, img_len=imgsize_arcmin, bin_k=512)
gmap_pwl_1024_k, gmap_pwl_1024_pk = fourier_power_spectrum(gmap_pwl_1024, img_len=imgsize_arcmin, bin_k=512)
/home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/fromnumeric.py:3440: RuntimeWarning: Mean of empty slice. return _methods._mean(a, axis=axis, dtype=dtype, /home/hergtl/.virtualenvs/py39env/lib/python3.9/site-packages/numpy/core/_methods.py:189: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(gmap_cmb_512_k, gmap_cmb_512_k**2/(2*pi)**2 * gmap_cmb_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(gmap_cmb_256_k, gmap_cmb_256_k**2/(2*pi)**2 * gmap_cmb_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the flat map $\\bf{generation}$ \nof the CMB power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(gmap_cmb_512_k, gmap_cmb_512_k**2/(2*pi)**2 * gmap_cmb_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(gmap_cmb_256_k, gmap_cmb_256_k**2/(2*pi)**2 * gmap_cmb_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_1024, label="nside=1024")
ax1.plot(gmap_pwl_512_k, gmap_pwl_512_k**2/(2*pi)**2 * gmap_pwl_512_pk, '-o', ms=2, c=c__512, label="nside=512")
ax1.plot(gmap_pwl_256_k, gmap_pwl_256_k**2/(2*pi)**2 * gmap_pwl_256_pk, '-o', ms=2, c=c__256, label="nside=256")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="Fourier power spectrum \nfrom the flat map $\\bf{generation}$ \nof the custom power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_1024)
ax2.loglog(gmap_pwl_512_k, gmap_pwl_512_k**2/(2*pi)**2 * gmap_pwl_512_pk, '-o', ms=2, c=c__512)
ax2.loglog(gmap_pwl_256_k, gmap_pwl_256_k**2/(2*pi)**2 * gmap_pwl_256_pk, '-o', ms=2, c=c__256)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
vmax = max(abs(fmap_cmb_1024[0].min()), abs(fmap_cmb_1024[0].max()))
plt.title("flat projection from full-sky CMB map\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_cmb_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
plt.figure()
vmax = max(abs(gmap_cmb_1024.min()), abs(gmap_cmb_1024.max()))
plt.title("flat map generation for CMB\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(gmap_cmb_1024, cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7fa20c10eb80>
vmax = max(abs(fmap_pwl_1024[0].min()), abs(fmap_pwl_1024[0].max()))
plt.title("flat projection from full-sky custom power-law map\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(fmap_pwl_1024[0], cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
plt.figure()
vmax = max(abs(gmap_pwl_1024.min()), abs(gmap_pwl_1024.max()))
plt.title("flat map generation for custom power-law\n %.f x %.f deg$^2$ \n %d x %d pixels"
% (imgsize_arcmin/60, imgsize_arcmin/60, pixsize_1024, pixsize_1024))
plt.imshow(gmap_pwl_1024, cmap=cmap, vmin=-vmax, vmax=vmax)
plt.colorbar()
<matplotlib.colorbar.Colorbar at 0x7fa20c4861c0>
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_proj, label="projected full-sky map")
ax1.plot(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_fgen, label="flat map generation")
ax1.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="CMB power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(fmap_cmb_1024_k, fmap_cmb_1024_k**2/(2*pi)**2 * fmap_cmb_1024_pk, '-o', ms=2, c=c_proj)
ax2.loglog(gmap_cmb_1024_k, gmap_cmb_1024_k**2/(2*pi)**2 * gmap_cmb_1024_pk, '-o', ms=2, c=c_fgen)
ax2.set_ylim(min(np.nanmin(fmap_cmb_1024_k[:20]**2/(2*pi)**2*fmap_cmb_1024_pk[:20])/10,
np.nanmin(gmap_cmb_1024_k[:20]**2/(2*pi)**2*gmap_cmb_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
fig = plt.figure(figsize=(fw, 1.5*fh))
ax1 = plt.subplot(211)
ax1.plot(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_proj, label="projected full-sky map")
ax1.plot(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_fgen, label="flat map generation")
ax1.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
ax1.legend(bbox_to_anchor=(1, 1), handlelength=0, labelcolor='linecolor', title="custom power-law power spectrum")
ax2 = plt.subplot(212)
ax2.loglog(fmap_pwl_1024_k, fmap_pwl_1024_k**2/(2*pi)**2 * fmap_pwl_1024_pk, '-o', ms=2, c=c_proj)
ax2.loglog(gmap_pwl_1024_k, gmap_pwl_1024_k**2/(2*pi)**2 * gmap_pwl_1024_pk, '-o', ms=2, c=c_fgen)
ax2.set_ylim(min(np.nanmin(fmap_pwl_1024_k[:20]**2/(2*pi)**2*fmap_pwl_1024_pk[:20])/10,
np.nanmin(gmap_pwl_1024_k[:20]**2/(2*pi)**2*gmap_pwl_1024_pk[:20])/10))
power_axes_k(ax1, ax2)
ax1 = plt.subplot(111)
ax1.loglog(fmap_cmb_1024_k, fmap_cmb_1024_pk, label="projected full-sky map")
ax1.loglog(gmap_cmb_1024_k, gmap_cmb_1024_pk, label="flat map generation")
ax1.set_xlabel("$k$")
ax1.set_ylabel("$P(k)$")
ax1.legend(handlelength=0, labelcolor='linecolor', title="CMB power spectrum");
ax1 = plt.subplot(111)
ax1.loglog(fmap_pwl_1024_k, fmap_pwl_1024_pk, label="projected full-sky map")
ax1.loglog(gmap_pwl_1024_k, gmap_pwl_1024_pk, label="flat map generation")
ax1.set_xlabel("$k$")
ax1.set_ylabel("$P(k)$")
ax1.legend(handlelength=0, labelcolor='linecolor', title="custom power-law power spectrum");